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A one-dimensional method for reconstructing the structure of prestellar and proto- 
stellar clouds is presented. The method is based on radiative transfer computations 
and a comparison of theoretical and observed intensity distributions at both mil- 
limeter and infrared wavelengths. The radiative transfer of dust emission is modeled 
Q H . for specified parameters of the density distribution, central star, and external back- 



ground, and the theoretical distribution of the dust temperature inside the cloud 
is determined. The intensity distributions at millimeter and IR wavelengths are 
computed and quantitatively compared with observational data. The best-fit model 
parameters are determined using a genetic minimization algorithm, which makes it 



possible to reveal the ranges of parameter degeneracy as well. The method is illus- 
trated by modeling the structure of the two infrared dark clouds IRDC-320. 27+029 
(P2) and IRDC-321. 73+005 (P2). The derived density and temperature distributions 
can be used to model the chemical structure and spectral maps in molecular lines. 

1. INTRODUCTION 

Massive stars are one of the most important components of our Galaxy and other galaxies, 
since they determine the energy balance in the interstellar medium and are the principal 
source of the heavy elements. However, the formation of massive stars is understood less 
well than the formation of low- mass stars (M < 10 M Q ). It is not yet known whether massive 
stars form like low-mass solar-type stars or stars with different masses form in different ways 
(see [1]). 

The general scenario of solar-type star formation has been developed with a fairly high 
degree of confidence. It starts with the formation of a dense collapsing core in a molecular 
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cloud. During the core collapse, a central protostellar object and circumstellar torus appear, 
which are later transformed into a young star and a gasBTj"dust disk, respectively. The 
durations of the various stages of this process and the details of the transitions between 
them remain open problems. However, objects at different stages of star formation have 
long been observed: prestellar (starless) cores, class and I protostars, classical T Tauri 
stars (class II protostars) and weak-lined T Tauri stars (class III protostars). 

No similar sequence of formation stages has been developed for massive stars, as a result 
of both observational and theoretical difficulties. All regions of massive star formation are 
far from us, and so require observations with high angular resolution (the distance of the 
closest, in Orion, is approximately 500 pc; typical distances are of the order of 2 kpc or 
more). The high gas densities in these regions lead to high extinctions, making optical and 
even IR observations difficult. As the energy of massive stars influences the surrounding 
gas starting from the earliest stages of their existence, regions of massive star formation 
have very complex structure, which requires the application of three-dimensional models in 
theoretical studies. 

The detection of massive prestellar cores — the counterparts of low-mass prestellar cores 
- could be an important step toward understanding the nature of massive star formation. 
We expect them to have the simplest structure, while still providing important information 
concerning the initial conditions for massive star formation. The most promising candidate 
objects are Infrared Dark Clouds (IRDCs). These clouds are seen in absorption against 
the Galactic IR background over wavelengths from several /im to several tens of /im; they 
were detected as a result of surveys with the ISO [2] and MSX [3] space telescopes. Their 
cores, i.e., the densest IRDC regions, can also be detected in emission at millimeter and 
submillimeter wavelengths. Some of these demonstrate signs of star formation [4], but the 
crucial test will be the detection of massive starless cores, i.e., without embedded compact 
sources and with collapse signature in their spectra. 

This makes searches for molecules whose spectral transitions could trace the kinematics 
of massive prestellar cores topical. Finding such molecules and transitions requires model- 
ing the chemical evolution and radiative transfer based on available information about the 
distributions of the gas and dust densities and temperatures in these objects. In the case 
of low-mass cores, information about their structure is obtained from data on either their 
emission in the submillimeter and millimeter [5], or the absorption of the background stellar 
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emission [6] in the optical and NIR. 

Spectra of infrared dark clouds can be studied in more detail. First, due to the presence 
of the infrared background, they can be observed in absorption over a much broader wave- 
length range than typical clouds in regions of low-mass star formation. Second, like other 
gasBTj"dust clouds, IRDC cores can be observed both in absorption and emission, making it 
possible to construct their spectral energy distributions (SEDs) in more detail, and, conse- 
quently, to reproduce their physical structure more reliably. However, this requires detailed 
numerical modeling and is a fairly resource-intensive task. 

In this paper, we present a method for studying the density and temperature distributions 
in prestellar cores, based on modeling the radiative transfer in them taking into account both 
their own radiation and the absorption of the background radiation at NIR to millimeter 
wavelengths. As an example, we use two dense IRDC cores from the sample of Vasyunina 
et al. [7], IRDC 320.23+0.32 and IRDC 321.71+0.07. The observational data at 1.2 mm 
described in [7] 1 are available for these cores, as well as Spitzer maps at wavelengths from 
3.5 to 70 /iin, obtained as a result of the GLIMPSE [8] and MIPSGAL [9] surveys. 

We selected the sources P2 in IRDC 320 and P2 in IRDC 321 for our study. Both these 
sources appear round in millimeter maps (Fig. 1), suggesting that a spherically symmetric 
approximation may be applicable. There is an important difference between the sources: 
70 /iin emission is associated with IRDS 321, whereas IRDS 320 is not detected in either 
emission or absorption at 70 /im. If we are dealing with massive prestellar or protostellar 
cores, the presence of emission at 70 /im may indicate a later evolutionary stage. For brevity, 
we further call these cores IRDC 320 and IRDC 321. Table 1 presents their coordinates, 
distances, and the masses and column densities derived from millimeter observations. 

The special feature of the method used here is that we reproduce the observational data 
at millimeter, near-IR, and mid-IR wavelengths in the framework of a single model. Similar 
modeling had been performed earlier, but using smaller numbers of wavebands. In particular, 
observations at 450 and 850 /im were used to determine the thermal structure of an IRDC 
core seen toward the giant molecular cloud W51 [10]. However, as we will show below, the 
use of only long wavelength data may be insufficient for this problem. 

1 Available at the Strasbourg astronomical Data Center (CDS) at http://cdsweb.u-strasbg.fr/cgi- 
bin/qcat?J/A+A/499/149 
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Figure 1. Intensity distributions in the IRDC 320 (left column) and IRDC 321 (right column) at 8 
/im (top row) and 70 /im (bottom row). The grey scale shows the IR emission and the contours show 

the emission at 1.2 mm. 

2. SEARCH FOR BEST-FIT CLOUD PARAMETERS 



The reconstruction of a source structure based on observational data is an inverse problem. 
Such problems are usually ill-defined and cannot be solved without some initial assumptions 
about the structure of the studied region. Usually, a specified model with several free 
parameters is used, which can be found by fitting the model to the observational data. 
When the number of free parameters is large, some optimization algorithm must be used to 
find the best fit. Such algorithms often implement a search for the minimum of a function of 
several variables, where some criterion for the agreement between the model and observations 
serves as this function and the relevant free parameters are the variables. Here, we searched 
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Table 1. 

Parameters of the modeled IRDC cores from [7]. 



Source a (J2000.0) 


5 (J2000.0) 


Distance, Mass (1.2 mnV 


I, Column density, 






kpc 


M 


10 22 cm" 2 


IRDC 320.23+0.32 (P2) 15 h 07 m 56.7 s 


-57° 54' 27" 


1.97 


50 


1.5 


IRDC 321.71+0.07 (P2) 15 h 18 m 26.7 s 


-57° 21' 56" 


2.14 


110 


3.2 



for the best-fit parameters applying the PIKAIA numerical code [11], which uses a genetic 
algorithm. This code yields both the localization of the functional minimum in the space of 
several parameters and the degeneracy of the parameters. Note that PIKAIA has already 
been used in some astrophysical studies (see., e.g.[12, 13]). 

2.1. Example of Applying the Method: A Two-Component Model 

We illustrate the method used to search for the best fit using the example of a simple two- 
component model of a protostellar cloud; various modifications of this procedure are fairly 
popular when analyzing observational data [14, 15]. We suppose that the cloud consists 
of two components, each with its own dust temperature T and dust surface density E. 
Physically these two components may be, for example, the cloud core and an envelope, but 
their physical interpretation does not matter from the computational point of view. For 
definiteness, we consider the first component to be located behind the second. The radiation 
intensities of the first and second components propagating toward the observer are then 
given by 

hly)= (l-e-^JS^TO (1) 
h{v) = e' K ^ h{v) + (1 - e- K - E2 ) B U (T 2 ), (2) 

where k v is the dust absorption coefficient per unit mass of dust and h^v) is the intensity 
received by the observer. These formulas neglect scattering of the light, and take into account 
only the intrinsic thermal emission of the medium. Thus, the radiation intensity in the two- 
component model depends on the four parameters Ti, T 2 , Ei and E 2 . For convenience, the 
surface densities Ei and E 2 can be converted into molecular-hydrogen column densities Ni 
and N 2 assuming the dust-to-gas mass ratio to be 0.01 [16]. The model will be compared 
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Figure 2. Results of searching for best- fit parameters using the two-component model for the 
protostellar cloud IRDC 321. The upper plot shows the history of the modelBT5™s convergence. 
The lower left plot demonstrates the localizations of the model parameters in T\-N\ space, and the 
lower right plot shows the same in T 2 -N 2 space. The black squares show models with x 2 < 1 and 

the dots show other models. 

with the SED observed over six wavelengths (1.2 mm, 70, 24, 8, 5.8, and 3.6 /im), using the 
observed intensities toward the assumed center of the cloud, determined as the peak of the 
1.2-mm emission. The chosen criterion for agreement between the model and observations 
takes the form 

1 N 

x 2 = ^£(M bs -ig'r d ) i°l (3) 
i=i 

where J° bs and 7™ od are the observed and theoretical radiation intensities for the ith frequency 
channel and ex, is the dispersion of the logarithm of the observed intensity. 

Figure 2 shows the results of the minimization based on the PIKAIA algorithm. The 
computation begins with arbitrary parameter values (only their limits are specified); there- 
fore, there is initially poor agreement with the observations, i.e., the x 2 values are very large 
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Figure 3. Comparison of the SED of the best- fit two-component model with observations of 
IRDC 321. The thin curves demonstrate the contributions of the first and second components, and 
the thick curve the total SED. The circles show the observed intensities. 

(Fig. 2, upper plot). After calculating about a thousand models, the algorithm finds a mini- 
mum with x 2 < 1; however, the search for the solution continues, since the algorithm checks 
whether the minimum is local or not. Further, there are no changes, even after computing 
10 000 models. The lower plots in Fig. 2 show the locations of the model parameters in the 
Ti-Ni and T 2 -N 2 spaces. The black squares represent the models with % 2 < 1, i.e., models 
that correspond to the minimum in the upper plot of Fig. 2. Figure 3 shows a comparison 
of the best-fit model with the observational data for IRDC 321. The model parameters 
are sharply localized and provide a good agreement between the theoretical SED and the 
observed points. 

Formally, the comparison of the model with the observations shows that the real cloud 
can be represented by a cold core with temperature Ti = 20 Pjb and gas column density Ni = 
2 x 10 22 cm -2 , surrounded by a warm envelope with T 2 = 400 Pjb and N 2 = 2 x 10 14 cm -2 . 
However, a comparably good agreement with the observed SED is given by a model with the 
opposite order of the components, i.e., with a hot core and cold envelope, and it is impossible 
to choose between these two models using this approach. 

Another shortcoming of the two-component model is its phenomenological character: we 
specify layer temperatures and column densities without considering how they could be 
realized physically, if they can be realized at all. Therefore, the results of such modeling 
should be interpreted with caution. For instance, the following questions arise in connection 
with the two-component model. First, if a cold core is surrounded by a warm envelope, how 
was the envelope warmed? Suppose that the cloud is warmed by some external radiation 
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field. If so, is it possible to choose parameters of this field that can reproduce the observed 
spectrum of the cloud? Second, if we suppose in a model with a hot core and cold envelope 
that the inner parts of the cloud are heated by some source, e.g., a protostar, is it possible 
to find protostar parameters that yield the necessary thermal structure of the cloud and 
reproduce the observed SED? 

We must bear in mind that this modeling took into account only the spectrum toward 
the center. However, we know both the SED toward the cloud center (or the integrated 
spectra) and the distribution of the radiation intensity, i.e., maps of the object at the different 
wavelengths, from the observations. In particular, the NIR intensity distributions for the 
IRDC cores show that the intensities decrease toward the cloud centers; i.e., these objects 
are seen in absorption, in contradiction with the initial assumptions of our two-component 
model. Thus, although such simplified models are formally able to describe the observed 
spectrum, they may not help us make progress in reconstructing the physical structures of 
protostellar clouds. 

Therefore, it seems logical to develop more physically feasible models of protostellar clouds 
and compare them with observations in more detail. A model of this type is described in 
the next subsection. 

2.2. Spatial Model of a Protostellar Cloud 

We suppose that the cloud is spherically symmetric, and its density distribution can be 
described, like those in low-mass prestellar cores [17], using the expression 



where n is the central number density of hydrogen nuclei, r is the radius of an inner 
region with approximately constant density n, and j3 is the index for the density decrease 
in the envelope. The outer radius of the cloud is taken to be 1 pc. In the radiative transfer 
modeling, we must also specify the inner radius of the cloud, which is always taken to be 



To take into account possible differences in the evolutionary status, e.g., the existence of 
a protostar, we assume that there is a source with radius i?* at the center of the cloud, which 




(4) 



50 AU. 
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emits as an ideal black body with temperature T*. The exact nature of this source is unim- 
portant: this may be the emission of either the protostar or accreting matter. The central 
source determines the temperature of the inner parts of the cloud. The cloud is irradiated 
from the outside by isotropic interstellar background emission with color temperature T^ g 
and dilution -Dbg, which determines the temperature in the outer parts of the cloud. 

In addition, we specify the isotropic IR background at the wavelengths used for compari- 
son with the observations. The intensity of this background is equal to the observed intensity 
at the edge of the cloud. At wavelengths shorter than 10 /im, the IR background is proba- 
bly due to the emission of small interstellar dust grains, polycyclic aromatic hydrocarbons 
(PAHs) [18]. Note that the intensity of the NIR background is higher than the intensity 
of blackbody emission with temperature Tb g and dilution D^ g , and significantly affects the 
cloudBTj™s appearance in the observations. 

The main agent determining the thermal structure of the cloud is dust, which absorbs, 
scatters, and reemits the continuum radiation. Therefore, we must solve for the continuum 
radiative transfer in order to find the dust temperature distribution inside the cloud and 
compute the distribution of the radiation intensity. 

We used the Accelerated Lambda Iteration (ALI) method for the radiative transfer model- 
ing, in which the mean radiation intensity is determined by integrating the radiative transfer 
equation along randomly chosen directions. The algorithm used is similar to that described 
in [19], but applies modifications necessary for thermal radiation. The temperature of the 
medium T is found from the equation of radiative equilibrium: 

oo oo 

J a v J v dv = J a u B u (T)du, (5) 



where a u is the absorption coefficient, J v = (47r) _1 J I u (n)dQ is the mean radiation intensity, 

4-7T 

I u is the spectral intensity of the radiation, and B v is the Planck function. In integrating 
the radiative transfer equation, 

(nV) I v = k v (S u - Q , (6) 

scattering is taken into account in the approximation of isotropic coherent scattering, when 
the source function S v takes the form 

b v = - . {() 

1%,, 



10 



In these equations, k v = a u + a u is the extinction coefficient and o v the scattering coefficient. 
The absorption and scattering coefficients as functions of frequency were computed according 
to the Mie theory for amorphous silicate grains using a code presented by D. Semenov 
(Max-Planck Institute for Astronomy, Heidelberg, Germany). We assumed that the size 
distribution of the dust grains is described by a power law, f(a) oc a~ 3 5 [20], with minimum 
and maximum grain radii of 0.001 and 10 /xm. 

When the temperature distribution is found, we can compute the theoretical distribution 
of the radiation intensity. For this, we used the temperature of the medium and the mean 
radiation intensity, which were determined when modeling the thermal structure of the cloud. 
To compare the theoretical and observational distributions of the radiation intensity, we 
convolved the theoretical distributions with the relevant telescope beams for each wavelength. 
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Figure 4. Results of radiative transfer modeling both with and without scattering for a 
representative protostellar cloud model. The left plot shows the temperature distribution over the 
cloud, and the right plot shows the intensity distribution at 24 fim. 

We will now demonstrate that we must include scattering in the radiative-transfer mod- 
eling. Figure 4 shows the model temperatures and intensity distributions at 24 /im obtained 
for a representative cloud with a central density no = 10 7 cm -3 and a central source tem- 
perature T* = 5000 Pjb, both with and without scattering. The temperature distribution in 
the cloud essentially does not depend on the scattering, but scattering significantly affects 
the distribution of the emergent radiation. Without scattering, the intensity distribution 
has deep minima, with the minimum depth determined by the column density toward the 
cloud center. After taking scattering into account, the decrease in the intensity toward the 
cloud center grows weaker, with the minimum intensity being determined by the ratio of the 
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absorption coefficient to the scattering coefficient. 

The method used to compute the temperature and intensity distributions was tested 
using a number of model problems. In particular, our computational results on the dust 
temperature agree with the solution obtained with the TRANSPHERE-1D numerical code, 
which was developed by C. P. Dullemond for modeling radiative transfer in dust envelopes 
taking into account the thermal emission of the dust 2 . 

3. MODELING RESULTS FOR IRDC 320 AND IRDC 321 

The method used to model the structures of protostellar clouds and the intensity distri- 
butions was applied to IRDC 320 and IRDC 321. We selected five variable model parameters 
(Table 2); three of these (n , r and /3) characterize the density distribution, and the other 
two (T* and -Dbg) the intrinsic and external radiation fields. The other parameters are fixed, 
and are also presented in Table 2. The inner radius of the cloud was specified in order to 
approximately take into account the absence of dust near the internal source (due to subli- 
mation). In principle, the dust sublimation radius depends on the source parameters, but 
we chose a fixed radius of 50 AU, which in our calculations corresponds to the upper limit 
for the size of the sublimation zones around hot stars. If there is no star, than the inner 
radius of the cloud essentially does not affect the distribution of emitted radiation. The outer 
radius of the cloud, 1 pc, is equal to the observed radii of the studied cores. The radius of 
the central source (star), 5 R Q [21], is chosen fairly arbitrarily, assuming that, within the 
explored range of stellar parameters, the heating will depend on the starBT)™s luminosity, 
and an incorrect choice of the starBTj™s radius can be compensated for through a suitable 
correction of its temperature. The temperature of the interstellar radiation field was taken 
to be 10 4 K. 

For each parameter combination, we modeled the radiative transfer, computed the inten- 
sity distributions, and quantitatively compared the computed and observed distributions. 
The search for the best-fit values of the variable parameters was done with the PIKAIA 
genetic algorithm. The models were compared with observations at 1.2 mm and 70, 24 and 
8 jim. As the clouds were assumed to be spherically symmetric, the radiation intensity 



2 http://www.mpia.de/homes/dullemon/radtrans 
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depends only on the angular distance from the cloud center. When constructing the one- 
dimensional observed intensity distribution, we took the center to coincide with the peak of 
the 1.2 mm emission. 

The agreement between the observed and model intensity distributions at each wavelength 
can be characterized using the standard \ 2 criterion. The generalized agreement parameter 
that must be minimized by the genetic algorithm is equal to the sum of the x 2 values for 
the individual wavelengths. We aim to determine the physical structure of the object using 
observations at both IR and millimeter wavelengths. The need to use both these wave bands 
is demonstrated by the ambiguity of the results obtained using only long-wavelength data. 
The modeling of IRDC 320 using only data at 1.2 mm formally yields the following best- 
fit parameters of the cloud: the central H2 density is 5.1 x 10 5 cm~ 3 , the plateau radius 
1.9 x 10 4 AU, the density-profile index 3.2, and the mass 580 M . However, this minimum 
is very flat. Figure 5 exhibits the model distribution in the T*-n(H.2) plane. We can see that 
the algorithm cannot distinguish between the models with and without an internal source. 
As an illustration, we selected two models denoted in Fig. 5 by the black squares. These 
are located at the edges of the range of allowed models: a less dense core with an internal 
source and a denser core without such a source. Figure 6 displays the distributions of the 
temperature and intensity for these two models. Even with the very different temperature 
distributions, the difference between the 1.2 mm intensity profiles is negligible. Thus, it 
is undesirable to use only millimeter wavelengths when determining the temperatures and 
densities of prestellar and protostellar cores. 

Further, we present our modeling results obtained using four wavelengths. The local- 
ization of the allowed parameters for both sources is shown in Fig. 7, which displays all 
models with \ 2 < 12 for IRDC 320 and x 2 < H f° r IRDC 321. The corresponding best- fit 
parameters for IRDC 320 and IRDC 321, are given in Table 2, together with the masses and 
hydrogen column densities toward the cloud centers. The temperature and column-density 
distributions for the best-fit models are shown in Fig. 8. 

In the case of IRDC 321, the allowed models are located in the region of higher densities, 
with a concentration of the density toward the cloud center (demonstrated by the high values 
of /3). IRDC 320 is characterized by flatter matter distributions and lower central densities. 
However, one must recognize that there remains some ambiguity in the solutions, even after 
fitting the model SEDs at four wavelengths. 
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Figure 5. Allowed models (\ 2 < 2.5), obtained as a result of fitting the theoretical intensity 
distribution for IRDC 320 to the observational data at 1.2 mm only. 
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Figure 6. Comparison of the radial temperature (left) and intensity (right) distributions for the 
two IRDC 320 models denoted by squares in Fig. 5. The dashed curve represents the model with 
central density 2 x 10 5 cm~ 3 and central source temperature 7300 K, and the solid curve the model 
with central density 1.2 x 10 6 cm~ 3 and without a central source. The spectra were fitted to the 

observational data at 1.2 mm only. 

The formal best fit for IRDC 321 (when the computation was terminated) has the central 
density 1.8 x 10 7 cm -3 , central plateau radius 5000 AU, (5 = 3.8, and central source tem- 
perature T* = 9300 K. However, there is another solution having nearly the same x 2 , with 
a higher central density of 1.8 x 10 8 cm -3 , central plateau radius of 2200 AU, and central 
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Figure 7. Localization of the allowed parameters for the IRDC 320 (circles) and IRDC 321 
(triangles) models in n - T* space (left) and R - (3 space (right) 
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Figure 8. Distributions of the density (left) and temperature (right) in the models for IRDC 320 

(solid) and IRDC 321 (dashed). 

source temperature of 6700 K. Thus, there remains some ambiguity for this source, which 
prevents us from distinguishing between models with denser, compact cores and colder stars 
and those with less dense, extended cores and hotter stars. Nevertheless, neither of the 
best-fit models admits the absence of a central source. In other words, the SED features 
of this source, especially the existence of the 70 /xm emission, cannot be explained if it is a 
starless source. 

The situation with IRDC 320 is slightly more definite. At the termination of the compu- 
tations, the best agreement with the observations was achieved for the model with central 
density 1.1 x 10 7 cm~ 3 , central plateau radius 4000 AU, f3 = 3.1, and central source tem- 
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Table 2. 

Parameters of the IRDC 320 and IRDC 321 models 

Parameter Notation IRDC 320 IRDC 321 



Adjustable parameters 

Central density of H 2 , cm" 3 n 1.1 x 10 7 1.8 x 10 7 

Plateau radius, AU r 4 x 10 3 5 x 10 3 

Index of density profile /3 3.1 3.8 

Star temperature, K % 7300 9300 

Background dilution D hg 1.3 x 10~ 13 6.7 x 10" 

Derived parameters 

Cloud mass, M M 170 230 

Column density, cm" 2 N 8.1 x 10 23 1.6 x 10 24 

Star luminosity, L Q L 60 160 

Fixed parameters 

Internal cavity radius, AU i?; n 50 

Cloud radius, pc R out 1 

Star radius, Rq R* 5 

Background temperature, K Tb g 10 4 



perature T* = 7300 K. Other models with similar x 2 values have approximately the same 
parameters; in contrast to the situation for IRDC 321, there is no alternative group of al- 
lowed models. Note that our fitting of the SED of this source over four wavelengths also 
indicates the existence of a central source. 

To test this conclusion, we determined the best-fit parameters for IRDC 320 excluding the 
heating by the central source, and obtained a worse agreement with the observed SED. In 
this case, we found that 70 jum is the crucial wavelength for determining whether the central 
source is present or not. The profiles of the source integrated intensity at other wavelengths 
can be reproduced equally well (or poorly) as in the model with the central source. 

Figure 9 compares the observed and theoretical intensity distributions for the best-fit 
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models. On the whole, a good agreement was achieved for both sources at all wavelengths. 
The models for IRDC 320 can explain the 1.2 mm emission, the absorption profiles at 8 
and 24 /xm, and the flat intensity distribution at 70 /im. The models for IRDC 321 yield 
emission at both 1.2 mm and 70 /im, whereas absorption profiles appear at both 8 and 
24 /im. The cloud masses are comparable, but the central density and column density is 
appreciably higher in the IRDC 321 model than the IRDC 320 model. The higher central 
density, column density, and stellar temperature in the IRDC 321 model compared to those 
in the IRDC 320 model lead to a higher 1.2-mm intensity and an increase in the emission at 
70 fim. The dashed curves in Fig. 9 show the intensity distributions calculated for clouds 
with the same parameters but without the central sources. The presence of a central star 
alters the intensities at 1.2 mm and 70 /xm, whereas the 8 /zm and 24 /zm intensities are not 
sensitive to the presence of this source. 

4. DISCUSSION 

We have described our method for studying the structure of dense molecular clouds based 
on modeling their SEDs at both millimeter and IR wavelengths. A number of studies recon- 
structing the density and estimating the masses of prestellar cores have used observations 
of dust emission only in the millimeter [5]. Since the intensity of this radiation in the op- 
tically thin approximation is proportional to the product of the dust column density and 
temperature, the dust temperature must be known to estimate the mass and density. 

Including shorter- wavelength (IR) data makes it possible to obtain a self-consistent re- 
construction of the density and temperature distributions. If the studied object is seen in 
absorption, i.e., the contribution of the dust thermal emission is negligible at these wave- 
lengths, the absorption intensity does not depend on the dust temperature, and is determined 
only by the dust surface density. On the other hand, if the cloud temperature is high enough 
to generate the intrinsic IR emission, the intensity of this radiation places strong constraints 
on the dust temperature. Supplementing these data with millimeter wavelength data signif- 
icantly narrows the range of allowed model parameters for the cloud. 

In other words, the millimeter-wavelength data must be supplemented with shorter- 
wavelength observations. However, modeling of the IR spectra encounters difficulties. First, 
the shorter the wavelength, the greater the importance of scattering. This does not seri- 
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ously affect the temperature distribution, but significantly alters the shape of the SED and 
the distribution of the integrated intensity (as was noted above). Including the effect of 
scattering appreciably increases the complexity of the models and makes calculations more 
resource-intensive. Further, the contribution of PAHs may be important at wavelengths of 
several /im. To take this into account correctly, we must make use of the parameters of both 
the external UV radiation and the PAH particles, which also increases the model complexity. 

In the current study, we included scattering in our modeling, but neglected possible PAH 
emission. Even in this form, the model makes us able to successfully reproduce a significant 
part of the observations at the four wavelength ranges used. Among shortcomings of the 
models, we note the poor fit quality at 24 /on: the central intensity minima in the models 
of both sources are deeper than the observed minima. The higher observed intensity may 
be related to scattering between the cloud and observer, which is not taken into account 
in our model. However, efficient scattering at 24 /im requires large dust grains, which are 
unlikely to exist in the interstellar medium. Another possible explanation is PAH emission 
in the cloud envelope (due to the influence of interstellar UV radiation). However, this 
explanation is also unlikely: the strongest PAH bands are in the 7-8 /im band, whereas 
our model successfully reproduces the intensity of the 8 /im emission. An important factor 
that could affect the spectrum of the outer parts of the cloud in the NIR is the stochastic 
heating of small dust grains [22]. Although including this effect should not significantly 
affect the internal structure of the cloud, it may be useful for determining the UV radiation 
intensity at the cloud periphery. Taking the foreground into account is also very important. 
In particular, strictly speaking, our conclusion that massive prestellar cores at 70 /im should 
be seen in absorption is valid for the cores considered only if all the background emission 
arises behind them. This assumption may be approximately correct for the relatively nearby 
clouds IRDC 320 and IRDC 321, but the contribution of foreground emission may be a critical 
component of the model in the more general case, and should be determined together with 
other parameters. 

Some of the computed models successfully reproduce the intensity profile at 24 /im, but 
not the distributions at other wavelengths. These are models with either extended, low- 
density cores (densities of the order of 10 5 cm" 3 ), or very compact cores with rarefied en- 
velopes. This may suggest a more complex morphology of the studied clouds, in particular, 
dumpiness and deviations from spherical symmetry. Another reason for the discrepancy be- 
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tween the model and the observations may be the failure of the assumption that the central 
source is compact. In reality, the cloud interiors may harbor a group of protostars, or may 
undergo volume heating as a result of gravitational collapse. We must also bear in mind the 
significant uncertainty in the distances of IRDCs, which affect the reconstruction of their 
geometric parameters. The method used here can be used to find the dust temperature 
distribution. However, modeling of the gas-phase chemical processes requires knowledge of 
the gas temperature rather than the dust temperature. The gas temperature can be directly 
determined from analysis of molecular line intensities; however, in this case, we encounter 
an ambiguity related to the complex chemical structures of the studied clouds and non-LTE 
conditions of the line formation. In particular, the temperatures obtained using ammonia 
lines are relevant only to those parts of the clouds where the ammonia molecules exist. Mod- 
eling of the dust temperature allows us to determine this quantity over the whole cloud. For 
the purpose of chemical modeling, we can assume that the gas and the dust temperatures 
are the same. This assumption is valid for the dense regions of molecular cloud cores [23], 
but fails at their periphery. Note, however, that the chemical composition at the periphery 
is essentially determined by photo-processes, whose rates do not depend on the temperature. 

In spite of the difficulties noted above, the results obtained demonstrate the effectiveness 
of our proposed method for reconstructing the physical structures of massive protostellar 
clumps. Further, we will use this method to model the chemical structure of these clumps 
and construct spectral maps, making it possible to better establish their evolutionary status. 
For this purpose, we must take into account the decoupling of the gas and dust temperatures 
at the periphery of the clumps, where gas heating due to the photoelectric effect becomes 
important. The low densities in these regions decreases the efficiency of gas cooling due to 
collisions with dust grains, thereby increasing the importance of gas cooling via molecular- 
line and atomic-line radiation, when the gas equilibrium temperature is higher than the dust 
temperature. 
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Figure 9. Comparison of the model and observed intensity distributions for IRDC 320 (left) and 
IRDC 321 (right). The solid curves show the results for the best-fit models, and the dashed curves 
show the results for the same models without the central sources. 
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